import pandas as pd
import numpy as np
# Visualisation libraries
## Text
from IPython.display import display, Markdown, Latex
## plotly
from plotly.offline import init_notebook_mode, iplot
import plotly.graph_objs as go
import plotly.offline as py
from plotly.subplots import make_subplots
import plotly.express as px
import warnings
warnings.filterwarnings('ignore')
In these methods, a smaller interval that contains a root is identified iteratively. As these methods do not require any derivatives in general, they are also known as no derivative methods. Some of examples for these methods are:
This method is used for estimatiion of the solution of the equation $f(x) = 0$ for the real variable $x$. Here, $f$ is required to be a continuous function defined on a closed interval $[a, b]$ and $f(a)$ and $f(b)$ have opposite signs. In this case, by the intermediate value theorem, $f$ must have at least one root in the interval $(a, b)$ (for further details see [1-2]).
def bisection(f, a, b, Nmax = 1000, TOL = 1e-4, Frame = True):
# endpoint values a, b,
# tolerance TOL,
# maximum iterations NMAX
if f(a)*f(b) >= 0:
print("Bisection method is inapplicable .")
return None
# let c_n be a point in (a_n, b_n)
Nmax=Nmax+1
an=np.zeros(Nmax, dtype=float)
bn=np.zeros(Nmax, dtype=float)
cn=np.zeros(Nmax, dtype=float)
fcn=np.zeros(Nmax, dtype=float)
# initial values
an[0]=a
bn[0]=b
for n in range(0,Nmax-1):
cn[n]=(an[n] + bn[n])/2
fcn[n]=f(cn[n])
if f(an[n])*fcn[n] < 0:
an[n+1]=an[n]
bn[n+1]=cn[n]
elif f(bn[n])*fcn[n] < 0:
an[n+1]=cn[n]
bn[n+1]=bn[n]
else:
print("Bisection method fails.")
return None
if (abs(fcn[n]) < TOL):
if Frame:
return pd.DataFrame({'an': an[:n+1], 'bn': bn[:n+1], 'cn': cn[:n+1], 'fcn': fcn[:n+1]})
else:
return an, bn, cn, fcn, n
n=n+1
if Frame:
return pd.DataFrame({'an': an[:n+1], 'bn': bn[:n+1], 'cn': cn[:n+1], 'fcn': fcn[:n+1]})
else:
return an, bn, cn, fcn, n
Example: Finding a root of the following polynomial. $$f(x)=x^{3}-x-2.$$
First, two numbers $a$ and $b$ have to be found such that $f(a)$ and $f(b)$ have opposite signs. Given
# defining function f using Lambdas
f = lambda x: x**3 - x - 2
For the above function,
a=1
display(Latex(r'$(a, f(a))$ = (%.2f, %.2f)' % (a,f(a))))
b=2
display(Latex(r'$(b, f(b))$ = (%.2f, %.2f)' % (b,f(b))))
Now since $f(a)<0$ and $f(b)>0$, we can implement the method. Let $N_{max}=20$ and the tolerance be $10^{-4}$
def Figs(data, Name, YL, CL = 'MediumBlue'):
fig = make_subplots(rows=1, cols=2, horizontal_spacing = 0.02, column_widths=[0.6, 0.4],
specs=[[{"type": "scatter"},{"type": "table"}]])
# Left
fig.add_trace(go.Scatter(x = data.index, y = data.fcn, mode='lines+markers', marker_color= CL), 1, 1)
fig.update_xaxes(range= [data.index[0]-0.08, data.index[-1]+0.08],
showgrid=True, gridwidth=1, gridcolor='Lightgray',
showline=True, linewidth=1, linecolor='Lightgray', mirror=True,
zeroline=True, zerolinewidth=1, zerolinecolor='Black', title_text = '$n$', row=1, col=1)
fig.update_yaxes(range= YL, showgrid=True, gridwidth=1, gridcolor='Lightgray',
showline=True, linewidth=1, linecolor='Lightgray', mirror=True,
zeroline=True, zerolinewidth=1, zerolinecolor='Black', title_text = '$f(c_n)$', row=1, col=1)
fig.update_layout(legend_orientation='v', plot_bgcolor= 'white', height= 600)
fig.update_traces(marker_line_color= 'Black', marker_line_width=0.5, opacity=1)
fig.update_layout(legend=dict(x=.02, y=.98, font=dict(size=12, color="black"), bordercolor="Blue", borderwidth=1))
fig.update_layout(title={'text': '<b>' + Name + '<b>', 'x':0.5, 'y':0.90, 'xanchor': 'center', 'yanchor': 'top'})
# Right
T = data.copy()
T.iloc[:,:-1] = data.iloc[:,:-1].applymap(lambda x: '%.6f' % x)
T.iloc[:,-1] = data.iloc[:,-1].map(lambda x: '%.2e' % x)
Temp = []
for i in T.columns:
Temp.append(T.loc[:,i].values)
Cols = [r'$a_{n}$', r'$b_{n}$', r'$c_{n}$', r'$f(c_{n})$']
fig.add_trace(go.Table(header=dict(values = Cols, line_color='darkslategray',
fill_color='DarkGreen', align=['center','center'],
font=dict(color='white', size=14), height=25),
cells=dict(values=Temp, line_color='darkslategray',
fill=dict(color= ['white','white','white', 'HoneyDew']),
align=['center', 'center'], font_size=12, height=20)), 1, 2)
if Name:
fig.update_layout(title={'text': '<b>' + Name + '<b>', 'x':0.5, 'y':0.90, 'xanchor': 'center', 'yanchor': 'top'})
fig.show()
data = bisection(f, a, b, Nmax = 20, TOL = 1e-4)
#
Figs(data, Name = 'Bisection method', YL = [-.5, 2])
Note that the difference between $x_n$ and a solution $x$ is bounded by
$$|c_{n}-c|\leq {\frac {|b-a|}{2^{n}}}.$$Therefore, The number $N_{max}$ of iterations needed to achieve a required tolerance $\epsilon$ is bounded by
$$N_{max}\leq \left\lceil \log _{2}\left({\frac {|b-a|}{\epsilon }}\right)\right\rceil ,$$def bisection_nmax(a, b, TOL): return int(np.ceil(np.log2(np.abs(b-a)/TOL)))
In our example, we need at least iterations.
bisection_nmax(a, b, TOL = 1e-4)
14
Thus, we can modify our algorithm and remove Nmax from it.
def bisection_mod(f, a, b, TOL = 1e-4, Frame = True):
# endpoint values a, b,
# tolerance TOL,
if f(a)*f(b) >= 0:
print("Bisection method is inapplicable .")
return None
# let c_n be a point in (a_n, b_n)
Nmax= int(np.ceil(np.log2(np.abs(b-a)/TOL)))+2
an=np.zeros(Nmax, dtype=float)
bn=np.zeros(Nmax, dtype=float)
cn=np.zeros(Nmax, dtype=float)
fcn=np.zeros(Nmax, dtype=float)
# initial values
an[0]=a
bn[0]=b
for n in range(0, Nmax-1):
cn[n]=(an[n] + bn[n])/2
fcn[n]=f(cn[n])
if f(an[n])*fcn[n] < 0:
an[n+1]=an[n]
bn[n+1]=cn[n]
elif f(bn[n])*fcn[n] < 0:
an[n+1]=cn[n]
bn[n+1]=bn[n]
else:
print("Bisection method fails.")
return None
if (abs(fcn[n]) < TOL):
if Frame:
return pd.DataFrame({'an': an[:n+1], 'bn': bn[:n+1], 'cn': cn[:n+1], 'fcn': fcn[:n+1]})
else:
return an, bn, cn, fcn, n
n=n+1
if Frame:
return pd.DataFrame({'an': an[:n+1], 'bn': bn[:n+1], 'cn': cn[:n+1], 'fcn': fcn[:n+1]})
else:
return an, bn, cn, fcn, n
Testing the above example,
data = bisection_mod(f, a, b, TOL = 1e-4)
#
Figs(data, Name = 'Bisection method', YL = [-.5, 2], CL = 'DarkGreen')
This method is a classic method that estimates a solution using the x-intercept of the line segment joining the endpoints of the function on the current interval. In this method, the root of the equation $f(x) = 0$ for the real variable $x$ is estimated when $f$ is a continuous function defined on a closed interval $[a, b]$ and where $f(a)$ and $f(b)$ have opposite signs. Moreover, $c_n$ is identified:
$$c_{n}={\frac {a_{n}f(b_{n})-b_{n}f(a_{n})}{f(b_{n})-f(a_{n})}}.$$def Regula_falsi(f, a, b, Nmax = 1000, TOL = 1e-4, Frame = True):
# endpoint values a, b,
# tolerance TOL,
# maximum iterations NMAX
if f(a)*f(b) >= 0:
print("Regula falsi method is inapplicable .")
return None
# let c_n be a point in (a_n, b_n)
Nmax=Nmax+1
an=np.zeros(Nmax, dtype=float)
bn=np.zeros(Nmax, dtype=float)
cn=np.zeros(Nmax, dtype=float)
fcn=np.zeros(Nmax, dtype=float)
# initial values
an[0]=a
bn[0]=b
for n in range(0,Nmax-1):
cn[n]= (an[n]*f(bn[n]) - bn[n]*f(an[n])) / (f(bn[n]) - f(an[n]))
fcn[n]=f(cn[n])
if f(an[n])*fcn[n] < 0:
an[n+1]=an[n]
bn[n+1]=cn[n]
elif f(bn[n])*fcn[n] < 0:
an[n+1]=cn[n]
bn[n+1]=bn[n]
else:
print("Bisection method fails.")
return None
if (abs(fcn[n]) < TOL):
if Frame:
return pd.DataFrame({'an': an[:n+1], 'bn': bn[:n+1], 'cn': cn[:n+1], 'fcn': fcn[:n+1]})
else:
return an, bn, cn, fcn, n
n=n+1
if Frame:
return pd.DataFrame({'an': an[:n+1], 'bn': bn[:n+1], 'cn': cn[:n+1], 'fcn': fcn[:n+1]})
else:
return an, bn, cn, fcn, n
Thus, for the above example, we have,
data = Regula_falsi(f, a, b, Nmax = 20, TOL = 1e-4)
#
Figs(data, Name = 'Regula falsi', YL = [-1, .2], CL = 'RoyalBlue')
There have been also some improved versions of this method. For example, in the Illinois algorithm, $c_n$ is identifed as follows.
$$c_{n}=\dfrac {{\frac {1}{2}}f(b_{n})a_{n}-f(a_{n})b_{n}}{{\frac {1}{2}}f(b_{n})-f(a_{n})}$$def Regula_falsi_Illinois_alg(f, a, b, Nmax = 1000, TOL = 1e-4, Frame = True):
# endpoint values a, b,
# tolerance TOL,
# maximum iterations NMAX
if f(a)*f(b) >= 0:
print("Regula falsi (Illinois algorithm) method is inapplicable .")
return None
# let c_n be a point in (a_n, b_n)
Nmax=Nmax+1
an=np.zeros(Nmax, dtype=float)
bn=np.zeros(Nmax, dtype=float)
cn=np.zeros(Nmax, dtype=float)
fcn=np.zeros(Nmax, dtype=float)
# initial values
an[0]=a
bn[0]=b
for n in range(0,Nmax-1):
cn[n]= (0.5*an[n]*f(bn[n]) - bn[n]*f(an[n])) / (0.5*f(bn[n]) - f(an[n]))
fcn[n]=f(cn[n])
if f(an[n])*fcn[n] < 0:
an[n+1]=an[n]
bn[n+1]=cn[n]
elif f(bn[n])*fcn[n] < 0:
an[n+1]=cn[n]
bn[n+1]=bn[n]
else:
print("Bisection method fails.")
return None
if (abs(fcn[n]) < TOL):
if Frame:
return pd.DataFrame({'an': an[:n+1], 'bn': bn[:n+1], 'cn': cn[:n+1], 'fcn': fcn[:n+1]})
else:
return an, bn, cn, fcn, n
n=n+1
if Frame:
return pd.DataFrame({'an': an[:n+1], 'bn': bn[:n+1], 'cn': cn[:n+1], 'fcn': fcn[:n+1]})
else:
return an, bn, cn, fcn, n
data = Regula_falsi_Illinois_alg(f, a, b, Nmax = 20, TOL = 1e-4)
#
Figs(data, Name = 'Regula falsi (Illinois algorithm)', YL = [-.15, .1], CL = 'Purple')
Interpolate Truncate and Project (ITP) method is an improved bisection method that can achive a superlinear convergence of the secant method. In worst-case, it matches the performance of the bisection method.
def ITP_Method(f, a, b, TOL = 1e-4, N0 = 1, K1 = 0.1, K2 = 2, Frame = True):
# endpoint values a, b,
# tolerance TOL,
# maximum iterations NMAX
Nmax = int(np.ceil(np.log2((b-a)/(2*TOL))) + N0)
# let c_n be a point in (a_n, b_n)
an=np.zeros(Nmax, dtype=float)
bn=np.zeros(Nmax, dtype=float)
cn=np.zeros(Nmax, dtype=float)
fcn=np.zeros(Nmax, dtype=float)
# initial values
an[0] = a
bn[0] = b
n = 0
while (bn[n]- an[n] > (2*TOL)):
mid = (an[n] + bn[n])/2
r = TOL * (2 ** (Nmax - n)) - (bn[n]- an[n])/2
xf = (an[n]*f(bn[n]) - bn[n]*f(an[n])) / (f(bn[n]) - f(an[n]))
sigma = np.sign(mid - xf)
delta = K1*((b-a)**K2)
#Truncation:
if delta <= np.abs(mid - xf):
xt = xf + delta * sigma
else:
xt = mid
# Projection
if np.abs(xt - mid) <= r:
cn[n] = xt
else:
cn[n] = mid - sigma * r
# Updating Interval:
fcn[n] = f(cn[n])
if fcn[n] > 0:
an[n+1]=an[n]
bn[n+1] = cn[n]
elif fcn[n] < 0:
an[n+1] = cn[n]
bn[n+1]=bn[n]
else:
an[n+1] = cn[n]
bn[n+1] = cn[n]
n += 1
if Frame:
return pd.DataFrame({'an': an[:n], 'bn': bn[:n], 'cn': cn[:n], 'fcn': fcn[:n]})
else:
return an, bn, cn, fcn, n
Furthermore, for the above example, we have,
data = ITP_Method(f, a, b, TOL = 5e-4)
#
Figs(data, Name = 'ITP Method', YL = [-.6, .6], CL = 'OrangeRed')